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AUTOMATIC ANALOG COMPUTER SCALING 
USING DIGITAL! OPTIMIZATION TECHNIQUES 


by 

John Celmer and Mary Rouland* 
Goddard Space Flight Center 


INTRODUCTION 

The Automatic Scaling Program is a digital program designed to eliminate the tedious, time- 
consuming process of manually scaling a linear or nonlinear system of differential equations for 
an analog computer. Optimum scaling performed by the program minimizes the deviation from 
unity of each amplifier gain. Since the analog's accuracy is limited to four significant figures, any 
scaling that results in excessively high or low gains, that is, with a large deviation from unity, re- 
duces precision even further. While it is important to recognize the program's time-saving fea- 
ture and optimization of the analog's limited precision, it is even more important to emphasize 
that true optimum scaling is extremely difficult if not impossible to perform manually. All the 
inputs required are simple numerical data; the outputs include documentation of all maximums, 
levels, and potentiometer readings and gains. 


GENERAL PRINCIPLES 

The Automatic Scaling Program is a digital program which optimizes the scaling of an arbi- 
trary number of differential equations, linear or nonlinear, for an analog simulation. The obvious 
advantage of having a digital computer perform the scaling is the time saved for the engineer who 
formerly performed this tedious task manually. The primary advantage, however, is that the pro- 
gram produces an optimum set of maximum values. 


The necessary input, designed for simplicity, requires only simple numerical data. The num- 
ber of equations, neq, the coefficient matrix, A(i, j), and an integer matrix m(I, J) which identifies 
the variable or variables of each term, suffice to completely describe a system to be scaled. The 
input also allows the engineer to specify some constant maximum values or to place minimum 
and/or maximum constraints on one or more maximum. If only neq, A(l, J), and m(I, J) are used 
as input, that is, if no maximums are to be held constant or constrained, the program assumes the 
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maximum value vector 


AMA(I) = 10.0 I = 1, 2, , NEQ 

as a starting point for the search for the optimum vector. 

Optimum scaling minimizes the deviation from unity of the gain at each amplifier. This gain, 
represented by a matrix, pg(I, J), describes the potentiometer gain of every term in all equations. 
Each element of the pg matrix 


PG(I, J) 


[A(I, J)j {ama[m(i, J)]} 
AMA(I) 


is calculated to minimize the merit value 


MERIT 


NEQ No. of terms 

[PG(I, J) - l] 2 

1=1 J=1 


L L 


where PG(I, J) is replaced by l/PG(l, j) if PG(I, J) < 1.0 and PG(I, j) is replaced by 1 if M(I, j) * i, 
in order that unchangeable first-order loop gains do not contribute to the merit value. 

The AMA vector associated with the minimum merit value of a set of equations, located by the 
program T s search routine running completely free without constant maximums or constraints as 
input, produces optimum scaling of the problem, that is, scaling that is subject only to problem 
coefficients. If some maximums are held constant or constrained, the scaling is optimum but it 
is subject to constraints and constant maximums as well as problem coefficients. 

The principal method used in obtaining the minimum merit value is the gradient or direction 
of steepest descent in conjunction with the Fibonacci search procedure— a one-dimensional, se- 
quential routine used to determine the optimum distance to move in the direction of steepest descent. 
In general, the slope or sensitivity of each maximum increases as its magnitude decreases, with 
very small maximums having very steep slopes. Therefore, if a problem has any maximums less 
than unity with steep slopes, the step size or distance moved in the steepest-descent direction is 
limited and the steepest-descent routine becomes ineffective in its convergence to the minimum 
merit value. The program follows a new search technique in which a modified steepest-descent 
movement is performed in three steps. With all AMA T s less than or equal to unity held constant, 
the first step allows maximums greater than unity to move in the modified steepest descent direc- 
tion until the merit value is minimized. The second step allows a similar movement for maximums 
greater than one-tenth but less than or equal to one, and the third step follows the same procedure 
for maximums less than or equal to one-tenth. 

The three-step modified steepest-descent method continues its normal sequence unless the 
auxiliary search routine is activated by the manual turning on of a predesignated sense switch on 
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the computer console. The auxiliary search is a combination sequential, linear interpolation 
routine in which the slope of each maximum less than one is individually minimized with respect 
to all other AMA’s held constant. While this would be too time-consuming to be a part of the normal 
search procedure, it is useful as a Tr kicking" device near the beginning of the search or whenever 
the three-step modified steepest descent may be moving very slowly. If the search is actually 
near the true global minimum, the rate of change of the merit value due to the auxiliary search 
will not be much improved over the normal routine and will be much more time-consuming. How- 
ever, if the modified steepest-descent search is slowed by a local minimum, the auxiliary search 
will "kick" it out of the region and allow the normal search to become effective again. Therefore, 
one of the two methods— the modified steepest-descent or the manually controlled auxiliary routine- 
will continue the search for the AMA vector that minimizes the merit value, until the precision or 
stepping condition specified by the program is reached. 

All data used as input in the program, including any optional data, are printed for future ref- 
erence and documentation. The merit value is printed at the conclusion of each search step in 
order that the programmer may follow the progress of the search and, with this information, de- 
cide when the auxiliary search might be useful. After optimum scaling is determined, documenta- 
tion of results is provided. This includes the equation number, the maximum value, the level, the 
potentiometer readings and gains of each equation, and the minimum merit value. An additional 
feature in which all levels are rounded to one significant figure is also included as output, since 
convenient levels may be preferred. The merit value resulting from the rounding of levels is also 
given, to simplify evaluation of the rounding effect. 


STRUCTURE 

Required Input 

The standard notation for systems of algebraic equations is 


X(l) 


Z> 


i. J) * X(J) 


(I, J = 1, 2 - --N) , 

where C(l,j) isanNXN coefficient matrix. For linear differential equations the standard nota- 
tion is 


X(I) 



i. J) * X(j) . 


(I, J = 1. 2 •••N) . 


3 


Since analog simulations involve both summations and integrations it is convenient to mix both 
algebraic and differential equations 


X(I) 


z> 


I, K) * X(K) 


X(J) 


E 


C(J, K) 


K 


* X(K) , 


(1 = 1- 

• • NA; J = NA + 1 • 

••N) 

(K -- 1 • 

••N) . 



The N x n matrix C frequently contains a profusion of zero entries (the number of nonzero entries 
in any row is usually less than five). It is thus more convenient to use an alternate notation in- 
volving an N x 4 coefficient matrix A in conjunction with an n x 4 integer matrix m 


X(I) = 2^ A( J ’ K ) * x[M(I f K)] , 

K 

x(j> = J2 A(J ’ K) * X K J ’ K) ] • 

K 

(I = 1 • • • NA; J = NA + 1 • • • N) , 

(K = 1 ••• 4) . 

This approach reduces storage requirements and simplifies input preparation. 

Assume the following set of differential equations: 

X(l) = - 0.5 [x( 1)] +0.09 [X( 2 )] - 2.0 [x(3) * X( 1)] , 

X( 2) = + 1.0 [ X( 1 ) ] - 2.0 [X( 2)] , 

X(3) = - 0.055 [X( 1) * X( 2) ] - 0.8 [X(3)] . 

The input required to scale this problem is neq (the number of equations), an A(l, J) matrix 
(coefficients of each term), and an M(i, j) matrix (the variable or variables of each term), where 
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I = 1, 2 , • • • , NEQ and j = 1, 2, — , number of terms of (I) . 



NEQ 

= 3 


A( I, J) Matrix 




1 

2 

3 

1 

-0. 5 

+0.09 

-2.0 

2 

+ 1.0 

-2.0 


3 

-0.055 

i 

o 

00 


M(I, J) Matrix 




1 

2 

3 

1 

1 

2 

300 1* 

2 

1 

2 


3 

1002* 

3 



Note: (*) indicates the form in which multipliers are entered as data. The format in which the 
program accepts these data will be included under ’’Operating Instructions— Required Input.” 

Optional Input 

The engineer may provide an initial set of maximums describing the point in space where the 
search will begin and/or he may specify that some of the maximums are to be held at constant val- 
ues. The maximum value of any equation not included in the optional input is set to 10.0 before the 
search is begin. If the system is nonlinear, the engineer may want to place minimum and/or maxi- 
mum constraints on the maximums of some equations. This is also possible in the optional input. 
Format for the four variations of optional input is found under ’'Operating Instructions -Optional Input.” 

Optimum Scaling 

Each amplifier or integrator used in the analog mechanization of the equations has a gain as- 
sociated with each input. The coefficients of a problem uniquely determine the total loop gain through 
any n th - order loop of the system where 1 < n < NEQ. While amplitude scaling cannot alter the total 
gain through any given n th -order loop, it can insure against any excessively high or low gains by 
even distribution. If all n th -order loops were independent of each other, calculation of the total 
loop gains of order n and the n th roots of these total gains would suffice to describe optimum dis- 
tribution of gains through each loop. Since the n th -order loops are not independent, an optimization 
technique far too complex to perform manually is employed. 

Program Notation 

The notation used for the vector whose n components are the maximum values of the system 
is AMA(i), where i = 1, 2, 3, — , neq. The elements of nter(I) where I = 1, 2, — , neq are the 
number of terms of each equation. 


5 



The gain at each amplifier is represented by the matrix PG(I, J); that is, the potentiometer - 
amplifier gain of each term in every equation. Each element of the PG matrix is calculated by the 
following equation: 


PG(I, J) 


[A(I, J)] {ama[m(I, J)]} 

~~ AMA(I) 


(1-a) 


If the M(I, j) of any term indicates that two variables v x and v 2 are to be multiplied, then the equa- 
tion is interpreted as 


PCX l, J) 


[A(I, J)] [ama(Vj) • ama(v 2 )] 
AMA(I) 


(1-b) 


Optimum scaling, defined previously as minimizing the deviation from unity of the gain at each 
amplifier, is accomplished by minimizing the merit value, 0, described by 


0 


NEQ NTER(I) 

[PG( I , J) - l] 2 

1=1 J=1 



( 2 ) 


where (a) PG(i, J) is replaced by 1.0 /pg(I, j) if PG(I, j) < 1.0 and (b) PG(I, J) is replaced by 1.0 
if M(i, j) = I. Part (a) above assures that low gains will not be ignored in the minimization rou- 
tine; for example, a PG gain of 1/4 contributes the same value to 0 as a PG gain of 4. Part (b) says 
that pot gains required in first-order feedback loops are not to be included in the value of 0. Since 


A( I , J) AMA[M(I, J)1 
PG(I,J) = AMA(Tr 


A( I , J) AMA(I) 
AMA(I) 


A(I, J) 


when M(i, J) = I, the pot gain is constant, not dependent on a maximum-value choice, and therefore 
not included in the calculation of 0 . 

Search Methods 

The principal search routine used in finding the minimum merit value is the gradient or 
steepest-descent method. The gradient is obtained from the derivative of the function 0 at some 
point in the space represented by the vector AMA 0 (I) , I = 1, NEQ. Since 0 is a multi-variable func- 
tion, (a function of NEQ variables), partial differentiation is used. The gradient line or direction of 
steepest descent is represented parametrically by 

<90 

AAMA o( X ) = <9A MA 0 (I ) * ( I - 1, 2, * * ■ , NEQ) , (3) 

where X is an arbitrary negative parameter. 
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The selection of an optimum \, which determines the distance to move on the steepest-descent 
line, results in movement to a new point in space, AMAj (I) , calculated from 


AMAj (I) = AMA 0 (I) + A AMA 0 (I) (I = 1 , 2, • • • , NEQ) . 


( 4 ) 


If the system of equations does not have product terms, the new merit value associated with 
the new point in space, AMAj (I) , is calculated from Equations 1-a, 2, 3, and 4 to obtain 


neq nter^i)J A(Ii j)|amA 0 [m(I, J)] + AAMA„ [ M ( It J)]} ) 

‘'AL " AMA 0 (I) +AAMA 0 (I) 1 


1=1 j=i ^ 


or simplified 


NEQ NTER 


*1 


L 


1=1 J=1 


V> [a(i, j> |amAj [m<i, j>]} _ 

AMAj (I) 


For a nonlinear system with multipliers, Equations 1-b, 2, 3, and 4 result in 

NEQ NTER 

*. ■ ll 


NEQ NTER | A( i , j ) [amA q ( V x ) + A AMA Q ( Vj )] [AMA 0 ( v 2 ) + A AMA 0 ( V 2 )] ^ 

. - 1 / 


i=i J=i 


AMA 0 (X) + A AMA 0 (I) 


(5-a) 


(5-b) 


(6 -a) 


or simplified 




NEQ NTER 


1=1 J=1 


A(I, J) [AMAj (Vj) ■ AMAj (v 2 )] 


AMA : (I) 



(6-b) 


where V x and V 2 are the variables to be multiplied. 

Since in both Equations 5 and 6 is a function of the independent variable \ in Equation 3 and 
has only one degree of freedom, a imi -dimensional search routine is used to find the value of X that 
minimizes <£ x . This routine consists of selecting trial values of k within the interval being searched, 
evaluating for these trial values, and reducing the interval. These three procedures take place 
according to a specified algorithm. 

Wilde (1964) discusses one- dimensional search procedures and includes a table showing re- 
duction ratios for various sequential search plans. The reduction ratios show the Fibonacci search 
to be the most efficient, followed closely by the Golden Section routine. The Fibonacci and Golden 
Section procedures are very similar once the first experiment (trial value) has been determined. 

The basic difference is that in the Fibonacci search the number of experiments is selected a priori 
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and enters into the calculation of the first experiment, while the starting point for the Golden Sec- 
tion is determined by a fixed ratio with the number of experiments dependent upon a defined error 
criterion. These techniques are effective in finding the minimum value of a function if it is not 
multimodal (more than one minimum) in the interval to be searched. If the function does have more 
than one minimum in the interval, it is possible that the search may not locate the minimum having 
the lowest value of the function. 

The search procedure used in the program to find the optimum \ is the Fibonacci technique 
using twenty experiments. The Fibonacci number associated with twenty experiments is 10,946; 
this says that after twenty experiments the original interval of uncertainty is reduced to less than 
0.0001 times its original length. 

The first step in the Fibonacci search is to establish a lower bound and an upper bound of the 
K interval. The lower bound, BDl, is set at 0 for all searches; the upper bound, BD2, is recalculated 
for each search. Since maximum values must be greater than zero and since X is always a nega- 
tive parameter, a unique upper bound, BD2, is determined before each search by the following: 


if 



(1 = 1,2,***, NEQ) , 


then 


BD2 * MINIMUM (Z r >0.0) (1 = 1,2,---, NEQ) . 


(V 


Also, Z r is negative only if the partial derivative, <?0/(9AMA o (I) , is negative. Re-examination of 
Equations 3 and 4 shows that maximums with negative partials cannot become negative and there- 
fore are eliminated in the calculation of BD2 in Equation 7. 

The next step is to select the first experiment, exp 1 , in the established interval. Wilde (pp. 29- 
30) gives the following formula for obtaining x 1? the ratio of the interval [BDl, exp 1] to the whole 
interval [BDl, BD2] : 


x 


i 


F "-i , ('I)"* 
F„ F„ 


(8) 


where n is the number of experiments, F n is the reduction ratio for n experiments from the table 
mentioned previously, and e is the minimum separation between any two experiments. Since the 
rate of change of 0 as a function of \ is very diverse, it would be difficult to choose a constant € 
appropriate for all searches. Computer experimentation showed that the formula, 
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was very close to the precision of Equation 8 with n = 20, and satisfactorily handled the varied 
slope field of the 4> function. 

The length of the interval to be searched is 

L = BD2 - BD1 ; (10) 

while the value for EXPERIMENT 1, the first trial value of k, is 

EXP 1 = BDl + L • Xj . (11) 


The search is continued by placing EXP 2 symmetric about the midpoint of the interval [bdi, BD2] 
with respect to EXP 1. Calculation of the <f> values with k = EXP 1 and \ = EXP 2 allows the ex- 
perimenter to discard a portion of the interval, either [BDI, EXP 2] or [EXP l, BD2] depending on 
which experiment has the minimum 4> value. In Figure 1 with the function <t> as shown, the segment 
[BDI, exp 2] would be eliminated from the interval since ^(EXP 2) > <£(EXP l) . Thus, the new BDI is 
EXP 2 and EXP 3 is placed symmetric about the midpoint of the new interval [BDI , BD2] with re- 
spect to EXP 1. After n experiments, the final 
selection of X is calculated from 

v _ BDI + BD2 

X - 2 • (12) 

With the optimum k known, Equations 3 and 4 
can be solved for the new point in the space, 

AMA 1 (i) , where I = 1, 2, • • • , NEQ. Following 
the same procedure, the vector Ama 2 (i) is 
found. This cycle continues until the new point 
is as close to the true optimum point as the 
precision of the program allows. 

In general, slope (sensitivity) of the merit 
value as a function of an AMA component varies 
inversely with the magnitude of the component. 

Calculation of X for any steepest-descent movement in which a given equation or equations have 
small maximum values and large positive partial derivatives shows that k must be extremely small 
in order to satisfy Equation 7, which assures that no maximums will become negative or zero. This 
obviously becomes a limiting factor on the rate of convergence to the optimum point. Even if all 
small maximums had negative partials, the size of k would still be limited owing to the steep slopes 
involved. 



PARAMETER 

Figure 1— Graph of merit value versus 
arbitrary negative parameter. 


Extremely small AMA's have very steep slopes. 



Since most linear and any nonlinear system of equations with a moderate range of coefficients 
necessitates presence of small maximums, a solution to the limitations of steepest descent is of- 
fered in the form of a three-step, modified steepest-descent movement: 

Step 1: a. AMA t s < 1.0 are held constant 

b. AMA t s > 1.0 move in modified steepest descent direction until 0 is minimized. 

Step 2: a. ama t s < 0.1 and ama t s > 1.0 are held constant 

b. 0.1 < AMA T s < 1.0 move as in step 1(b). 

Step 3: a. ama t s >0.1 are held constant 

b. AMA T s <0.1 move as in step 1(b). 

If at any time in the search, the 20-step Fibonacci routine does not provide enough precision 
for the merit value to be steadily decreased, the program automatically increases the number of 
steps of the X search. This feature is desirable, particularly when the neighborhood of the optimum 
point is reached. 

Although the modified steepest descent is fast and sufficiently handles most problems, an aux- 
iliary search was programmed which can be activated and deactivated by sense switch 5, turned on 
and off manually on the computer console. The auxiliary search is a combination sequential, 
linear-interpolation routine, in which the partial derivative or slope of each maximum less than one 
is individually minimized with respect to all other maximums held constant. While in most cases 
this is too time-consuming to be a part of the normal search, it is useful near the beginning of a 
search or at any time that the modified gradient method is moving slowly, perhaps indicating ap- 
proach to a local minimum. The auxiliary search provides a boost so that the modified gradient 
can again be effective. 

Special Systems 

A LCC (linear constant coefficient) system of equations allows a great deal of flexibility in the 
optimization of scaling. A suggested approach for automatic scaling is to use only the required in- 
put as discussed under "Structure— Required Input," letting the program supply the initial amA vector. 

Since all LCC systems have an infinite number of ama vectors associated with the optimum- 
merit test, the resulting ama vector may then be multiplied by any scalar and still maintain the 
same minimized merit value and corresponding potentiometer settings. This property of a linear 
system enables the engineer to keep peak voltages near the ±100-volt level simply by adjusting the 
scalar quantity. 

A linear system with variable coefficients, while not permitted the freedom of a LCC system, 
does have an infinite number of AMA vectors that have the same minimum 0 value. However, since 
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these maximum vectors are not proportional, as they are in the LCC case, the engineer must rely 
on the search routine to find different vectors with the same minimum <t>. 

Nonlinear systems have only one maximum- value (AMA) vector having the minimum merit value. 
This vector will be located by the search if no constraints or constant maximums are specified as 
input. If constraints are placed on some maximums or if any maximums are determined and used 
as input, scaling will be optimum but subject to the input constraints. If limits and nonlinearities 
other than multipliers are in the system, it is suggested that the engineer examine the nonlinearities 
and place constraints or constant maximums where they are necessary or beneficial. 


Output 

The output provides documentation of (1) all data used as input to the program, (2) the results 
of the optimum scaling— maximum values, levels, potentiometer readings and gains (printed se- 
quentially by equation), and (3) the minimized 4> value. 

Obviously, the levels chosen by the program, although optimum for the analog simulation, are 
not convenient rounded numbers. Since in most cases the engineer prefers scaling that will allow 
him to quickly change from scaled to problem units, perhaps even mentally, the output includes an 
additional feature in which all levels are rounded to one significant figure, barring any which had 
been specified as constant in the input. Revised documentation of (2) and (3) listed above is printed. 
The increased 4> value provides simple evaluation of the rounding effect. 

A data card is punched for each equation, giving the equation number and the rounded maximum 
value in the format accepted as input to the program. This deck is invaluable to the programmer 
for any future rescaling that may be required by additional constraints or changes or simply for 
new documentation. 

All input prepared for the Automatic Scaling Program (NEQ, A MATRIX, M MATRIX) and 
the AMA deck (punched output of the program discussed in the preceding paragraph) also is in the 
proper format for the FAST Program, a digital program which provides static and dynamic solu- 
tions as a checkout for the analog. 


OPERATING INSTRUCTIONS 
Required Input 

Assuming the three -equation problem with the equations, coefficients, and variables as dis- 
cussed under T 'Structure— Required Input,” the following cards are required. Data card 1 contains 
NEQ entered as an integer, right- justified in columns 1-5. Data cards 2-4 contain the A(i, J) ma- 
trix, one card for each equation. In card 2, for instance, columns 1-5 contain the equation number, 
right-justified, integer form. The coefficient A(l, 1) is entered in floating point starting in 
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/ 3 



/ 7 

- .5, .09, -2.0 

COL 1-5 

7-72 

COL 1-5 

7-72 


Card 1 


Card 2 


/ 2 


^0, -2.0 



r 

'.055, -.8 

COL 1-5 

7-72 

COL 1-5 


7-72 


Card 3 


Card 4 


column 6 and terminated by a comma. The column to the immediate right of the comma contains 
the sign or first digit of A(l, 2). A(i, 2) is similarly terminated by a comma. This process 
continues until all coefficients of the equation have been entered. Each coefficient including signs 
and decimal point must not exceed 15 columns. 



Card 5 


Data cards 5-7 contain the M(I, J) matrix, 
one equation per card. The equation number in 
interger form is right- justified in columns 1-10 
and is followed by the M(l, 1), M(I, 2), M(I, 3) 
terms , all i nteger, right- justified in columns 
11-20, 21-30, 31-40. 

Card 5 indicates the format for multipliers. 
The integer 3001 right-justified in the 10 col- 
umns allowed for the m( 1, 3) shows that the vari- 
able of term (1, 3) is the product, X 3 * . 



Card 6 


Card 7 


12 


I 




Optional Input 

The optional input is entered as two different data sets described by: 

Set 1: Contains a starting set of maximums and indicates which, if any, of these values are 

to be held constant. 

Set 2: Specifies all minimum and maximum constraints placed on any AMA ! s. 

A blank card must be at the end of each data set to indicate termination of the set unless values 
are given for every equation. Even if a data set is omitted entirely, the blank card is still needed. 

Assuming the same problem used as an example in the required input section, suppose that 
the programmer wants the maximums of Equations 1 and 2 to have the values 12.0 and 5.6 at the 
beginning of the search. Suppose also that the 5.6 value is to be held constant throughout the search. 
These values are found on cards 8 and 9. Columns 1-5 contain the equation number in integer form 
and right-justified. The maximum appears in floating point starting in column 6, not to exceed 
column 21. Any maximums to be held constant must have a minus sign preceding the equation 
number. 

Since Equation 3 is not supplied with a starting or constant maximum, the program sets its 
maximum equal to 10.0 before the search is begun. Card 10 must be a blank card to show the end 
of optional input, Set 1. 

Suppose that the AMA (1) should not exceed 100.0 and the ama( 3) should be greater than 5.0 but 
less than 50.0. These constraints are found on cards 11 and 12. The equation number is entered 
as an integer, right- justified in columns 1-5. The constraints are entered as floating-point num- 
bers; first the minimum which starts in column 6, next a comma, and finally the maximum. If 
there is no maximum constraint, the comma is not necessary. 


/ lM 2.0 


COL 1-5 7-72 


/ -2U.6 


COL 1-5 7-72 


Card 8 


Card 9 


/ 


it 

1 0 100 . 


COL 1-5 


7-72 


Card 11 




t 

35] 


50. 


COL 1-5 7-72 


Card 12 


13 


Card 13, a blank card indicating the end of optional input Set 2, is the final card of the input. 


Sense-Switch Operation 

Table 1 shows the programming controlled by sense switches 1-5. The switches marked TT * Tf 
are particularly valuable. Under normal conditions all switches are turned off. 


Table 1 

Programming Controlled by Sense Switches 1-5. 


s.s. 

Number 

’’ON 1 ' Operation 

1 

Prints documentation after each search. 

2 

Prints Fibonacci steps and Auxiliary Search. 

3 

Resets program for more than one set of data. 


Designed to allow the programmer to choose his own precision; 

"ON" state terminates the program, provides complete documentation. 

*5 

Activates auxiliary search. 


FUTURE DEVELOPMENTS 

Future developments will focus on additional flexibility and ease in programming nonlinearities. 

Division of variables will be allowed with standard input form similar to the present form of 
multiplier input. 

Beside permitting constant maximums and constraints as now programmed, future development 
will include the option of stating constant ratios to be maintained between two or more maximum 
values. This will be helpful in eliminating potentiometers and handling built-in gains in switching 
networks and other equipment. 

If mathematical time- scaling is beneficial, this must, at present, be precalculated and reflected 
in the coefficient matrix. In the future the engineer will be given the option to experiment with dif- 
ferent time-scaling factors by simply providing the factors as input. The program will evaluate 
the merit of time-scaling by comparing all minimized <t> values. 


Goddard Space Flight Center 
National Aeronautics and Space Administration 
Greenbelt, Maryland, August 1969 

604-31-03-12-51 


REFERENCES 

Wilde, D. J., ’’Optimum Seeking Methods," Englewood Cliffs, N. J.: Prentice Hall, Inc., 1964, 
pp. 10-52. 
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Appendix A 


Automatic Scaling Example 

The sample problem displayed in this section was run using the Fortran program listed in 
Appendix B developed specifically for the SDS 9300 digital computer. The sample program shows 
scaling of the set of three differential equations previously discussed in both sections dealing with 
required input. 


AUIO^mT lu SCALifM^ tXAMPLEi (NO OPTIONAL INPUT) 


>\ '*1 A 1 K I x 


1 

2 

6 


A • O v J 'J U J J u 
• jbb j 


•■jVUUUvUU 

“ 2 • u U U u u ^ u 0 
” • J'-'CvOu 


'2 • uu^-’UvOv 

.'JL 

■ v o v JUt Ul 1 


• UuUUwL-U'J 

• UJJUt'U'JU 


■•1 ,1 A I rV I A 


A i 

2 1 

3 1ov2 


2 ^Clwl 
2 

3 


■1 A X i *| u rl ^ 


0 o ,m b 1 K A i f,, 1 ^ 
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automat 1 c 

EG 
+ + + 

SCALING EXAMPLE 
AM A ( I > 

♦++++++++++++++ 

DELTA( I ) 
+++++++++++++♦+ 

levelii) 
•»■ + + + + + + + + + + ■»■ + + + 

0A 1 £ 
TERM1 

+++++++++♦+++++ 

is apr 1909 page 

TERM2 

+++++++++++++++ 

00U2 

TERM3 

+ + ■*■ + + + + + + + + 4- + + + 

TLRM4 

+ 4 + + + + + 4 4-4 + + ■»*♦ 

1 

10-000000 

• 00 

1 0* OuOuOO 

-50000o 

•09000U 

*0. 000000 


2 

lO.UuOOOO 

• 00 

IO.UUOUOO 

1 .OOOOOu 

2 * OOOOOU 



3 

IO.UuOuOO 

.00 

lu.uuouoo 

. 550U0U 

•800000 




MERIT 

TEST 


463- 9 ^ 3 ^ 393^9 

NO. 

1 

MERIT = 

47*6914760942 

N 6 . 

4 

MERIT = 

31 .8816575337 

NO. 

7 

mlk i t - 

2u. 1429973378 

NO. 

tu 

MERIT = 

I7.16t6107626 

NO. 

1 3 

MERIT = 

15.1862255049 

NO. 

16 

MERIT - 

14.6474871758 

NO. 

19 

mLRIT = 

14.3417070574 

NO . 

22 

MERIT = 

14.27^2215434 

NO. 

25 

MERIT = 

14.2284847689 

NO. 

28 

MERIT = 

14.2226775511 

NO. 

31 

merit = 

14.2155651171 

NO. 

34 

MERIT = 

14. 1997667058 

NO. 

37 

MERIT = 

14.0170847617 

NO. 

38 

MERIT = 

14.0161128830 

NO. 

40 

MERIT = 

12.8622146225 

NO. 

41 

MERIT = 

12.8622130664 

NO. 

43 

MERIT =* 

12.75w896738i 

NO. 

44 

MERIT = 

12.7477292105 

NO. 

46 

MERIT = 

12.7426480996 

NO. 

47 

MERIT = 

12.7420924535 

NO . 

49 

MERIT = 

12.741 1870731 

NO. 

5u 

MERIT = 

12.7410791704 

NO. 

52 

merit = 

12.74^8452642 

NO. 

53 

MERIT * 

12.74o8163005 

NO. 

55 

MERIT = 

12. 74^7442527 

NO. 

56 

MERIT - 

12. 74u7 348893 

NO. 

58 

MERIT = 

12.74^7102898 

NO. 

59 

MERIT * 

12.74^7^72995 

NO. 

61 

MERIT = 

12.7 4w69877 9 1 

NO. 

62 

MERIT = 

12-74 06975928 

NO. 

64 

MERIT = 

1 2 . 74^694 1 927 

NO. 

65 

MERIT = 

I2.74g6937588 

NO. 

67 

MERIT = 

12.74 W 6924722 

NO. 

68 

MERIT = 

12.74^6923282 


AUTOMATIC 

SCALING 

;• EXAMPLE 


NO * 

70 

MERIT * 

12. 74u69 1 907 0 

NO* 

71 

merit * 

12.74o6918531 

NO . 

73 

MERIT * 

1 2 • 7 4(j69 1 697 q 

NO . 

74 

NErIt = 

12.74U6916693 

NO. 

76 

MERIT = 

12.74W6915765 

NO. 

77 

MERIT = 

12.74Q6915749 

NO. 

79 

MERIT * 

12.7406915453 

NO. 

80 

MERIT = 

12.7406915699 

NO. 

80 

MERIT = 

12.74o69l5446 

NO. 

82 

MERIT * 

12.7406916012 

NO. 

82 

MERIT * 

12.7406915399 

NO. 

83 

MERIT = 

I2.74o6915650 

NO. 

83 

MERIT = 

12.74o6915393 

NO. 

85 

MERIT = 

12.7 40691 6 1 i 7 

NO. 

85 

MERIT = 

12.7406915374 

NO. 

86 

MERIT = 

12.7406915658 

NO. 

86 

MERIT ■ 

12:7406915372 

NO . 

88 

MERIT = 

12.7 4w6 9 16197 

NO. 

88 

MERIT = 

12.7 4u69 1 5364 

NO. 

89 

MERIT = 

12.7406915669 

NO. 

89 

MERIT * 

12.7406915363 

NO. 

89 

MERIT = 

12.7406915363 

NO i 

89 

merit = 

12.7406915363 

NO . 

91 

MERIT = 

12.74^6916267 

NO. 

91 

MERIT * 

12. 74 06915362 

NO . 

92 

MERIT « 

12.7406915647 

NO . 

92 

MERIT a 

12.7406915359 

NO . 

94 

MERIT = 

I2.74u69i6283 

NO. 

94 

MERIT = 

12. 74^6915358 

NO. 

95 

MERIT = 

12.7 4o6915684 

NO • 

95 

MERIT = 

12.7 4u69 1 5358 

NO. 

95 

MERIT = 

12.7406915358 

NO . 

95 

MERIT = 

12.7406915358 

NO . 

97 

MERIT = 

12.7406916315 

NO. 

97 

MERIT = 

12.7 4o69 1 5358 

NO. 

97 

MERIT = 

12.7406915358 

NO. 

97 

MERIT = 

12.74o6915358 

NO . 

98 

merit = 

l2.74o6915684 

NO . 

98 

MERIT = 

12.74o69l5358 

NO . 

98 

MERIT = 

12.74w6915358 

NO . 

98 

MERIT = 

12.7 406^ 1 5358 

NO . 

100 

MERIT = 

12.7406916315 

NO . 

100 

MERIT = 

12.7406915358 

NO . 

100 

MERIT = 

12.74U6915358 

NO. 

100 

MERIT = 

12.7406915358 


AJ fOMAT IC 
EU 
+ + + 

SCALING EXAMPLE 
AMA ( n 

+ + + + + ■*• + -M- + + + + + + + + 

DEL i A ( 1 ) 

LEVEL ( I ) 
++++++++++++++ 

1 

3.0354 74 

1.81 

32.943779 

2 

1 U * 692903 

1.53 

9.180290 

3 

.953589 

i 

e 

c 

lu4. 869206 


MERIT TEST = ]2. 7 4 069l53 6 

MINIMuM MER I T y ALjjE 


AUTOMATIC 

EU 

+ + + 

SCALING EXAMPLE 
AMa ( i ) 

delta ( i j 

LEVEL ( I ) 
+■ + + + + + + + + + + + + + 

1 

3*333333 

i-®i 

3 U . OuOUOO 

2 

11.111111 

1.53 

9 • OuQOUO 

3 

1-000000 

-.00 

luo-ouoooo 


MERIT TEST = l2.9b4334 7 Q5 

ALL LEVELS ARE ROUNDED 


*SfOP* QOUOIOGU 



2*037037 


800000 



Appendix B 


Fortran Program 


c 


INTEGER NEQ.NTER ( 150 >»M( 150.4 >»C( 150 > 

DOUBLE PRECISION A ( 1 50. 4 ) * AM A ( 05 1 50 > » CT ( 1 50. 2 > * DELT A ( 0 : 1 50 > » 

PG ( 150*4 I.SAMA! 150 >»SDELTA( 150 >»XAMA( 150 >,XDELTA 
( 150> * XLEVl 150>,YAMA<150>.YDELTA<150> 

DOUBLE PRECISION ABS. SD 1 * BD2, BD3. BM AX, CD£i-T A. D, D 1 * D2, D I FF, EXP 1 * 
EXP2»FACT,GK,MERIT»SGK,SlN,SIS,SMERIT,SUBT,SUM» 
SUMP»V»X.X1PHI,X2PHI,Y»Z,ZZ 




c . 

C AUTOMATIC ANALOG SCALING 

C 

**************************************************** '2 ******************* 

C SEC l - READ IN NEQ, COEFFICIENTS, AND M MATRIX 


ASS!X>=DAB31X> 

READ ( 105, 10) NEQ 
10 FORMAT ( I5,4(F16.8> ) 

OUTPUT ( 108 ) NEQ,".", 'A MATRIX'," 

DO 30 1=1, NEQ 

READ ( 105, 10) C(I). (Al I»J)»J=1,4> 

HR I TE ( 108, 10> C( I ), < A l I,J),J = 1»4) 

IF(C( I I.NE.I ) GO TO 1370 
NTER ( I )=1 
DO 20 J= 1 » 4 

IF(A(I,J).EG.O.) SO TO 30 
NTER ( I ) *NTER ( I >+l 
20 CONTINUE 
30 CONTINUE 

OUTPUT ( 108)"»"»*M MAIRIX'," 

DO 50 1=1, NEQ 

RE AD ( 105, 40 ) C(I).<M!l.J).J*l.NTER<I)-l) 

WRITE! 108, 40> C< I ). <M< I . J ) » J=1 . NTER ( I )-l ) 

40 FORMAT (5110) 

I F ( C ( 1 ) .NE. 1 ) GO TO 1370 
50 CONTINUE 

C SEC 2 - RESET CONDITIONS FOR NEW AMA SET AND CONSTRAINTS 

60 SENSE LIGHT 0,3,15 

NUM*MUM=NAC=NOD*LAST=MAST*0 
BD3* 1 » 

AMA (0 ) =1 ,0 
SMER1T=.9M 10. **20) 

CDELTA-O.O 
70 DO 80 1=1, NEQ 

SDELT A ( I >=0.0 
AMA! I ) *12345. 

C ( I >*+l 
DO 80 J* 1, 2 
80 CT ( I . J ) =0. 


C SEC 3 - READ IN INITIAL MAXIMUM VALUES 
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I II 


OUTPUT ( 108 )"»".' MAXI MUMS'," 
oe 100 1 = 1. NEQ 

READ ( 105* 90) N» V 
IF (N.EQ.O) SO TO 110 
WRITE1108. 90) N.V 
90 FORMAT! 15. F16. 8) 

L=IA9S(N> 

AM A ( L ) = V 
C!L)=N 

100 CONTINUE 

110 DO 120 1=1. NEQ 

IF (AMA( I ) .NE. 12345. )S0 TO 120 
AMA ( I )=10.0 
120 CONTINUE 

SEC 4 - READ IN CONSTRAINTS 

OUTPUT ( 108 ) ".".'CONSIRAINTS'," 

DO 140 1=1. NEQ 

RE AD ( 105. 130) N. <CT (N.J ) » J = 1 » 2 > 

IF (N.EQ.O) SO TO 150 

WRITE! 108. 130) N,(CT(N,J),J=i,2) 

130 FORMAT! I5.2CF16.8) ) 

140 CONTINUE 

SEC 5 - CALCULATION OF MERIT VALUE 

150 DO 170 1=1. NEQ 

DO 170 J=1.NTER! I >-l 
LL = M ( I » J ) / 1000 
NN=M0D(M( I, J), 1000) 

160 PS! I. J ) = ( AMA(LL)*AMA(NN) *A< I. J) ) / AMA ( I ) 

PG< I.J)=ABS!PG( I, J) ) 

170 CONTINUE 
MER I T=0 • 

DO 200 1=1. NEQ 

DO 200 J = l» NTER ! I >-l 
IF(M(I»J) .EQ. I ) GO TO 200 
IF(PG( I.J) .GE.1.0) SO 10 180 
Y = (l.O/PG! I.J) J-1.0 
GO TO 190 
180 Y=PG< I.J >-1.0 
190 MER I r = MER I T + ( Y**2) 

200 CONTINUE 

SEC 6 - PRINT AMAS. LEVELS. POT GAINS. MERIT VALUE 


1 F ( SENSE LIGHT 3) 260. 220 
CONTINUE 

IF (SENSE SWITCH 2) 230» 250 
WRITE(108. 240) LY. MERIT 
FORMAT <5X.$LY = $, I 3. $. MER I T = $. F25 . 9 ) 
IF! SENSE LIGHT 10) 1080. 110U 


L 


I 



260 IF ! SENSE LIGHT 15) 270» 370 
270 WRITEU08# 280) 

280 FORMAT! 1H1,5X.$EQS»11X.$AMA( I ) $# 10X* SDELT A ( 1 ) S» 10X. SLEVEL < I )$»13X, 
* ST£RM1$*13X*STERM2$, ^.STERNS*. 12X.STERM4S ) 

WRITEU08* 290) 

290 FORMAT !5X, 3 !S+$)»2X, ( 7 1 15 ( S + S > * 3X ) ) ) 

00 340 1 = 1* NEQ 

XLEVU )*100./AMA(I i 
0UTPUT1108) * * 

IF! <MAST.EQ.l).OR.!LASI.EQ.l) ) GO TO 300 
GO TO 340 

300 1F<<MAST.£Q.1).AND.<CU).LT.0>> GO TO 330 
310 WRITE) 106# 320) I.AMAU) 

320 FORMAT! 15.F16. 8) 

GO TO 340 
330 IT = -I 

WRITE! 106* 320) IT.AMAi I ) 

340 WRITE! 108* 350> I . AMA < I ) • SDELT A ! I ) , XLEV 1 1 ) * < PG ! I , J ) , J«1 . NTER 1 I ) - 1 ) 
350 FORMAT ( 2/* 5X» !3*2X»F15*6*F18»2*3X» < 5 < F 15 . 6* 3X ) ) ) 

WRITE! 108» 360) MERIT 
360 FORMAT !5/*5X»SMERIT TEST =$*F24.9> 

IFILAST.EQ.l) GO TO 1450 
I F < M AST .EQ. 1 > GO TO 1400 
OUTPUT ! 108 ) ' ' 

SENSE LIGHT 2 

370 l F ( SENSE LIGHT 2» 400* 380 
380 WRITE! 108» 390) NOD.MERIT 
390 F0RMAT!5X,$N0.$. I4.5X.SMERIT = S.F21.10) 
************************************************************************ 
C SEC 7 - COMPARE NEW MERIT WITH MINIMUM MERIT* STORE MINIMUM MERIT 

************************************************************************ 
400 IF (SENSE LIGHT 14) 1390* 410 

410 IF! (SMER1T-MERIT) .GT. ( 10. * * < - 10 ) ) ) GO TO 460 
IF1NAC.NE.2) GO TO 430 
DO 420 1 = 1* NEQ 
420 AMA! i ) = Y AMA I I ) 

NAOO 
GO TO 570 
430 NUM=NUM+1 

1F1NUM.LT. 4) GO TO 1050 
DO 440 1 = 1* NEQ 
440 AMA! I )=YAMA< I ) 

MUM=MUM+ 1 
NUM*0 
GO TO 510 

450 SENSE LIGHT 3*14*15.16 
GO TO 1250 
460 SMER I T=MER I T 

DO 470 I =1 »NEQ 
Y AM A < I ) *AM A ! I ) 

YDELTA1I )=S DELTA! I ) 

470 CONTINUE 
NUM»MUM=0 

IF! ! MOD! NOD* 3) .EQ.2) .AND. (NAC.EQ.2) )NOD=NOD+l 
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SEC 8 - SENSE SWITCH T3 TERMINATE PROGRAM 


IP ( SENSE SWITCH 4) 480* 490 
480 SENSE LIGHT 3*15 
MAST*1 
GO TO ISO 


C SEC 9 - SENSE SWITCH FOR INTERPOLATION OF AMAS LESS THAN ONE 


490 I F ( SENSE SWITCH 5) 500» 510 
500 NAC*l 

GO TO 520 
510 NAC“0 

IFIMUM.GT.3) GO TO 450 


C SEC 10 - CALCULATION OF DELTAS 


520 NOD*N0 D* 1 

530 1F!M0D!N0D.3).EQ.1> GO TO 570 
IF1M0DIN0D.3) .EQ.O) GO TO 550 
DO 540 I *1. NEQ 

IF((AMA(I).GT.0.1).AND.(AMA(1).LE.1.0>) go TO 570 
540 CONTINUE 
GO TO 520 

550 DO 560 1*1. NEQ 

I F ( AMA< I ) .LE.0.1 ) GO TO 570 
560 CONTINUE 
GO TO 520 

570 DO 560 1*1»NEQ 
DELTA! I >*0.0 
DO 580 J*1»NTER! I >-l 
NN*MOD(M( I.J). 1000 > 

LL*M( I.JI/1000 

IF(PG(I,J).LT.1.0> X*(2.*( (AMA( I > / < AM A < LL ) * AM A l NN ) * A < I , J > ) >-l. > >• 

* ( l./< AMA(LL>*AMA(NN)*A( 1. J) > > 

IF(PG( I, J > ,GE. 1.0> X=( 2.*! < ( AM A (LL >*AMA!NN>*A( I, J ) ) /AMA ( I ) ) -l .0 > > * 

* (-1.*! (AMA(LL)*AMA(NN)*A( I.J) >/ < AMA! 1 )**2) ) ) 

580 DELTA! I )*DELTA! I >+X 

DO 590 1*1. NEQ 

DO 590 J*t*NTER ( I )-l 
LL*M ( I , J ) 71000 
NN*MOD!M( I . J > . 1000 ) 

IF(PG(I,J).LT.1.0> Z*<2.*( ! AMA! 1 > / < AMA < LL > *AMA < NN ) * A < I.J) > )-l.) > 

* * ! ( - 1 . *AMA! I ) )/< A( l.J)*AMA!NN)*!AMA!LL>**2> > > 

IF IPG! I. J) .GE. 1.0) Z* ( 2. » l ( ! AMA ! LL ) * AMA < NN > * A ! I . J > ) /AMA ! I)) - 1 .0 ) > 

* * ( (A! I » J ) *AMA 1 NN ) ) /AMA I I ) > 

590 DELTA!LL>*DELTA<LL)+Z 

DO 600 1*1. NEQ 
DO 600 J*1»NTER ( I )-l 
LL*M! I. JI/1000 
NN=MOD ( M < I.J). 1000) 

IF!PG( I.J) .LT.1.0) ZZ*<2.* ( ( AMA! I )/( AMA ( LL ) * AMA ! NN ) * A ! I.J) ) )-l. ) ) 

* *!<-!. *AMA! I ))/!*( I.J )*AMA!LL)*< AMA (NN)*«2) ) ) 
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IF(PGd»J> .GE.1.0) ZZ=<2.*< < < AMA<LL)*AMA<NN)*Ad.J) ) / AM A ( I >>-!.»» 
* *< <AU*J>*AMA(LL) )/AMA(l) ) 

600 DELTA(NN)-DELTA(NN)+ZZ 


C SEC II - INTERPOLATION ROUTINE ACTIVATED dY SENSE SWITCH 5 


IF(NAC.E3.0> GO TO 970 
1F<M0D(N0D»3>.EQ.1) GO TO 970 
610 I F ( SENSE LIGHT 11) 730* 620 
620 I F ( SENSE LIGHT 19) 800* 630 
630 N I T-0 

REPEAT 900* WHILE ( < BM AX < AMA. DELTA, NEQ. C ) .G I . < . 01 » ) . AND . C N 1 T .LE. 20 > > 

DO 890 L = 1 * NEQ 

MA0=MAR=O 

S I S=S I N-. 01 

MM-MMM-0 

FACT-2.0 

SENSE LIGHT 0 

IF(C(L).LT.O) GO TO 890 

IF( AMA(L) .GT.l. ) GO TO 890 

IF ( ABS ( DELTA (L) ) .LE. ( . 01 ) ) GO TO 890 

X AMA ( L ) -AMA ( L ) 

XDELTAIL)-DELTA(L) 

640 AMA ( L ) -XAMA ( L ) 

I F ( SENSE LIGHT 12) 680* 650 
650 AMA(L)-AMA(L)+(AMA(L»*SIS) 

I F ( SENSE SWITCH 2) 660* 670 
660 OUTPUT ( 108 ) 'AMA + ,01'*AMA(L) 

670 GO TO 720 

680 AMA(L)«AMA(L)-( AMA(L)*SIN) 

690 1 F ( SENSE SWITCH 2) 700* 710 
700 OUTPUT ( 108 ) 'AMA - .Ol'.AMAU) 

710 SENSE LIGHT 17 
720 SENSE LIGHT 11 
GO TO 570 

730 I F ( SENSE SWITCH 2) 740* 750 
740 OUTPUT ( 108 IXDELTA (LI* DELTA (L I 
750 1F( (XDELTA(L)*DELTA(L) T.GT.O* 1 GO TO 820 
SUM=ABS(XDELTA(L) ) +ABS < DELT A < L I ) 

DIFF = ABS(XAMA(L)-AMA(L) ) 

SUMP= A8S(XDELTA(L) ) /SUM 
I F ( SENSE SWITCH 2) 760* 770 
760 OUTPUTdOS) XAMA(L),SUMP,DIFF 

770 IF(XAMAIL) .GT.AMA(L) ) AM A ( L ) *X AM A < L ) - < SUMP* D I FF ) 

IFtXAMA(L) .LT.AMAIL) ) AM A ( L ) -X AMA ( L > + ( SUMP* U I FF ) 

SENSE LIGHT 23 
1 F ( SENSE SWITCH 2) 780* 790 
780 OUTPUT1108) AMA ( L ) »L 
790 SENSE LIGHT 19 
GO TO 570 
800 MMM-MMM+ 1 

I F ( SENSE LIGHT 23) 810* 89C 

810 IF( ( ABSIDELTA(L) ) .LE. ( .1 ) ).0R. (MMM.NE.l ) ) Go TO 890 
IF(AMA(L) .GT.l.Q) GO TO 890 
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SENSE LIGHT 0 

mm=mad=mar=o 

SIS=SIN=.Ol 

FACT=2.0 

X AMA ( L ) »AMA ( L ) 

XDELTA<L)=DELTA<L) 

GO TO 650 

i IF(AMA(L).GE.CT(L#1) ) GO TO 630 
AMA < L ) =CT ( L» 1 ) 

C(L)=-L 
CT(L»1 )»0. 

OUTPUT(108) 'MINIMUM CONSTRAINT IN ITERATION' 
GO TO 790 

IF(AMA(L). GE.(. 00001 > > GO TO 840 

AMA(L)=l.E-5 

GO TO 790 

IF< AMA(L) .GE.l.Ol GO TO 790 
I F ( SENSE LIGHT 17) 860* 850 

IF( ABS(DELTA<L> ) . GE . ABS < XDELT A < L > ) ) GO TO 870 
MAR=MAR+ 1 

I F ( l MAR. EG. I ) .OR. (MAR.tG.2I ) SIS=SIS*10. 

IF (MAR.GE.3) SIS=SIS*2. 

GO TO 640 

IF<ABS<DELTA<D ) ,LT.ABS< XDELTA(L) > ) GO TO 880 
MM = MM+ 1 

IF ( MM . NE • 2 ) GO TO 880 
AMA(L)*=XAMA<L> 

GO TO 790 
SENSE LIGHT 12 
MAD=MAD+ 1 

IF t M AD . EQ • 2 ) SIN*SIN*10. 

IF (MAD.LE.2) GO TO 640 

IFIMAD.GE.5) FACT=1.0+‘ ( F ACT- 1 . 0 > /2 . 0 > 

S I N = S I N*F ACT 
GO TO 640 
CONTINUE 
N I T = N I T + 1 

I F ( SENSE SWITCH 2) 910* 920 
OUTPUT<108) '91 CONTINUE' 

SENSE LIGHT 3 
N AC*2 

I F ( SENSE SWITCH 1) 940* 950 
SENSE LIGHT 15 
DO 960 IM»NEQ 
SDELT A ( I ) = DELT A ( I ) 

GO TO 150 

SEC 12 - STORE AMAS ANU DELTAS 


970 DO 980 I=1#NEG 

SDELT A ( I ) “DELT A ( I ) 
S AM A ( I )=AMA< I > 

960 CONTINUE 
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SEC 13 - ESTABLISH UPPER BOUND OF LAMBDA 


990 MAX=0 

IF (M0D(N0D.3) • EQ . 2 > BD3M.E-5 
IF (M0D(N0D.3) . EQ ♦ 0 ) BD3M.E-8 
DO 1020 LXM.NEO 
lF(SDELTA(LX).LE.O.) GO TO 1020 
IF ( C ( LX ) .LT . 0 > GO TO 1020 

IF< <M0D(N0D.3> .EQ.l > .AND. (SAMA(LX) .LE.1.0* > GO TO 1020 
IF ( (MOD (NOD. 3) ,EQ.2» .AND. ( ( SAM A ( LX ) . GT . 1 . 0 > . OR . ( SAM A ( LX l.LE.O.l > > > 
* , GO TO 1020 

IF( <M0D(N0D»3) .EG.O) .AND. (SAMA(LX) .GT. .1 ) )G0 TO 1020 
1000 MAX*MAX + 1 

GK*ABS ( AM A (LX ) /SDELT A ( LX > > 

IF(MAX.NE.l) GO TO 1010 

SGK*UK 

GO TO 1020 

1010 IF(SGK.GT.GK) SGK*GK 
1020 CONTINUE 

IF(MAX.EO.O) BD2*BD3 
IF ( MAX.NE.O > BD2*SGK 
BDl s O. 

BD3*BD2 

l F ( SENSE SWITCH 2) 1030. 1050 
1030 WRITEdOS* 1040 > BD 1 . 8D2 
1040 F8RMAT(5X,$BD1*S.F20.9»$# BD2=$.F2Q.9> 


C SEC 14 - FIBONACCI SEARCH FOR OPTIMUM LAMBDA 


1050 DO 1190 L Y* 1 * 1 9 

IF(LY.NE.l) GO TO 1060 
D I FF*BD2"BD 1 
D2*10946. 717711. 

Dl*D2*DIFF 
EXP1*BD1+D1 
EXP2«BD2-D1 
SENSE LIGHT8 
SENSE LIGHT 9 

1060 1 F ( SENSE LIGHT 8> 1070» 1090 
1070 CDELT A*EXP1 

SENSE LIGHT 10 
GO TO 1250 
1080 X 1PH l =MER I T 

I F ( SENSE LIGHT 9) 1090* 1110 
1090 CDELTA-EXP2 
GO TO 1250 
1100 X2PHI=MERIT 

1110 IF( X2PHI .GT.X1PHI) GO 10 1150 
BD2*EXP1 
EXP1-EXP2 
SUBT«BD2-EXP1 
EXP2*BD1+SUBT 
X 1 PH I =X2PH I 

I F ( SENSE SWITCH 2) 1120. 1140 


1120 WRITE! 108* 1130) EXP 1 * LXP2* BDl * BD2 

1130 F0RMATi5X,$CASElS*Fi8.9*F18.9,$*BDl=$»Fl5.9.S,BD2* 5 S»F15.9) 

1140 I F ( L Y . EG . 19 ) GO TO 1200 
GO TO 1190 
1150 8D1«EXP2 
EXP2-EXP1 
SUBT-EXP2-BD1 
EXP1*BD2-SUBT 
X2PH1=X1PH1 
SENSE LIGHT S 

IF ( SENSE SWITCH 2) 1160* 1180 
1160 WRITE! 103* 1170> EXP 1* t XP2* BDl. BD2 

1170 F0RMAT(5X,$CASE2$»Fl6.9,F18.9»$*BDl=$»Fl5.9,$,BD2 ! *$»F15.9) 

1180 IF1LY.EG.19) GO TO 1200 
1190 CONTINUE 
1200 CDELTA* ! BD2+BD1 )/2.0 
SENSE LIGHT 3 

IF ! SENSE SWITCH U 1210* 1220 
1210 SENSE LIGHT 15 
1220 CONTINUE 

IF ! SENSE SWITCH 2) 1230* 1250 
1230 WRITE! 108* 1240) CDELTA 
1240 FORMAT <5X»$FINAL CDELT A=$, F2U. 9 > 

•••I********************************************************************* 

C SEC 15 - CALCULATION OF AMAS WITH OPTIMUM LAMBDA 

****************«******o************************************************ 

1250 DO 1340 I - 1 »NEQ 

1 F ( SENSE LIGHT 16) 1320* 1260 
1260 IF!C! I ) .LT.O) GO TO 1340 

IF1!M0DIN0D*3).EG.1).AND.!SAMA( I).LE.1.0»> GO TO 1340 
IF! 1 MOD! NOD* 3) .EQ.2) .AND. ! (SAMA! I ) . GT . 1 . 0 ) . OR . ( S AMA ( I ) .LE.0.1 > > > 
• GO TO 1340 

IF! !M0D!N0D*3).EQ.0).AND.!SAMA( I).GT..1))G0 TO 1340 

DELTA! I )*SD£LTA1 I )*!CDtLTA> 

1270 AM A < 1 ) *SAM A ! I ) -DELTA ( I ) 

IF! <CT!l,l>.EQ.0.>.AND.<CT(I»2>.EQ.0.)> GO TO 1340 
IF!(CT(I*1).NE.0.).AND.(CT(I*2).NE.0.)) GO TO 1290 
IFIC1 II.D.EG.O.) GO TO 1280 
IF! AMA! I ) .Gt.CT! I»l ) >G0 TO 1340 
GO TO 1300 

1280 IF! AMA! I ) .LE.CT! 1*2) > GO TO 1340 
GO TO 1310 

1290 IF! ! AMA! I ) .GE.CT! 1* 1 ) ) .AND. ! AMA! I ). LE.CT! 1 *2) > > GO TO 1340 
IF! AMA! I ) .GT.CT! 1*2) ) GO TO 1310 
1300 C ! I ) *- I 

AMA! I ) *CT 1 1 * 1 » 

CTlI.D-O. 

SENSE LIGHT 4 
GO TO 1340 
1310 C ( I > *-I 

AMA! I )»CT< 1*2) 

CT! I»2)=0. 

SENSE LIGHT 4 
GO TO 1340 


